Cluster algorithm for non-additive hard-core mixtures 
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In this paper, we present a cluster algorithm for the numerical simulations of non- additive hard- 
core mixtures. This algorithm allows one to simulate and equilibrate systems with a number of 
particles two orders of magnitude larger than previous simulations. The phase separation for sym- 
metric binary mixtures is studied for different non-additvities as well as for the Widom-Rowlinson 
model (B. Widom and J. S. Rowlinson, J. Chem. Phys. 52, 1670 (1970)) in two and three dimen- 
sions. The critical densities are determined from finite size scaling. The critical exponents for all 
the non-additivities are consistent with the Ising universality class. 

PACS numbers: 61.20. Ja, 64.60.Hr, 05.10.Ln 



I. INTRODUCTION 

The Widom-Rowlinson (WR) model attracted a lot 
of attention as a prototype for the liquid-vapor phase 
separation. This simple model is composed of a two- 
component system where likewise particles do not inter- 
act whereas unlike particles interact through a hard-core 
potential. It was shown to present a phase separation at 
high density and a critical point that belongs to the Ising 
universality class 0, • This mixture exhibits a liquid- 
liquid critical point (with large composition fluctuations) 
in contrast to pure fluids that experience a liquid-vapor 
critical point (with large density fluctuations). However, 
the phase transitions are related 3 and believed to de- 
pend on the same universality class [f| . 

A straightforward generalization of the WR model is 
the non-additive hard-core (NAHC) mixtures where the 
likewise particles also experience a hard-core interaction. 
Two particles i and j present a minimal distance of ap- 
proach ijxy with X and Y representing respectively the 
A or B component from which belong particles i and 
j. Non-additive mixtures are characterized by <jab — 
(c j 4 j 4+o'b_b)(1+A)/2 with A ^ 0. Additive mixtures cor- 
respond to A = 0. For a negative non-additivity A < 0, 
particles tend to form hetero-coordinations [|| 0. For 
a positive non-additivity A > 0, an entropically driven 
phase separation occurs between two phases at suffi- 
ciently high density due to the extra repulsion between 
unlike particles Q. The two phases are chemically dif- 
ferent, one is rich in A particles whereas the other is rich 
in B particles. This phase separation occurs even for the 
symmetric mixtures where a aa — ®bb- The critical den- 
sity as well as the universality class have been determined 
from different numerical simulations 0, 0, fToL [III IT^| . 
The special case of oaa = o~bb = corresponding to the 
WR model has also been studied 0,0. Additive mix- 
tures with a strong asymmetry a a a <C <jbb are also of 
interest and present another kind of entropically driven 
phase separation transition predicted for the first time 
by Biben and Hansen but only recently observed by 
numerical simulations |l4l Il5j . 

Most recent simulations of the NAHC mixtures 0, 0, 



Il2| used a Monte-Carlo algorithm within the semigrand 
canonical ensemble 0|. In addition to the simple 
moves of particles as simple Monte-Carlo steps, some 
steps consist in changing the nature (or component) of 
a particle when the hard-core interactions permit this 
modification. The resulting ensemble corresponds to a 
fixed total number of particles but a variable composi- 
tion (or fixed difference of chemical potentials between 
particles of different components). A detailed descrip- 
tion of the algorithm in the case of binary mixtures with 
squared-well interactions may be found in the paper of 
de Miguel et al. ^3- Working in the semigrand canon- 
ical ensemble gives access to the coexistence curve of 
the model. However, the algorithm is limited to rather 
small system sizes. The largest simulations concerned 
16384 particles [Til but most of them were limited to a 
few thousands |9l IfOl IT^ . Furthermore, only small non- 
additivities A < 1 have been simulated [ill E^ . The 
main reason for this limitation is the following. When 
a Monte-Carlo step corresponding to a change of com- 
ponents is tried, at high non-additivity there is a large 
probability of overlap with at least one of the neighbors. 
This results in the rejection of the change of components 
and a dramatic slowing down of the equilibration of the 
numerical simulations. 

In the present paper we consider a cluster algorithm 
proposed by Dress and Krauth 18] for hard-core mix- 
tures which proved useful for the detection of the phase 
separation in additive asymmetric mixtures [LH \l9L l20j| 
and for the analysis of two dimensional polydisperse 
hard-core mixtures [2l| . The advantage of this algorithm 
is to allow one to equilibrate systems with up to 10 6 par- 
ticles (two order of magnitude larger than previous sim- 
ulations). Also, the slowing down of the equilibration 
for large non-additivity is completely avoided with this 
cluster algorithm. This allows us to analyze the limit 
of infinite non-additivity (A — > oo) where the NAHC 
mixtures converge to the WR model. Furthermore, the 
coexistence curve is accessible in contrary to the previous 
cluster algorithm used for the WR model 

The paper is organized as follows: in section II, we 
present the model of non-additive hard-core mixtures and 
we describe the cluster algorithm used for the numerical 
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simulations. We present the numerical results for the 
two and three dimensional mixtures in section III. From 
finite size scaling analysis, we extract the critical densities 
and the critical exponents. Finally, we conclude by a 
discussion of the results in section IV. 



II. NON-ADDITIVE HARD-CORE MIXTURES 
A. Description of the model 

We consider non-additive hard-core (NAHC) mixtures 
in two and three dimensions. The system is made of two 
components A and B with respectively Na and N B parti- 
cles. The particles experience hard-core interactions. No 
overlap is possible if the distance d between the center 
of the particles i and j is lower than o~xy where X and 
Y represent the components of particles i and j respec- 
tively. For convenience, the particles are placed in two 
identical boxes of equal volume V. Periodic boundary 
conditions are assumed on each of the boxes. From the 
following description of the cluster algorithm, the reason 
for considering two boxes will become obvious. Notice 
that the consideration of two boxes is of common use for 
the determination of the coexistence curve in semigrand 
canonical ensemble simulations. 

In the following, we restrict ourselves to symmetric 
mixtures with a aa = ®bb = & however we allow for a 
general positive non-additivity A = ctab/c — 1- F)ue to 
the hard-core interactions, the temperature plays a triv- 
ial role and the phase diagram is determined only by the 
number density p = pa + Pb = Na/2V + N B /2V and the 
composition xa — ^ ~ xb — Pa/ (pa + Pb) of the system. 
In the particular case of the symmetric NAHC mixtures, 
the critical point (p c ,x c ) is determined in composition 
(x c = 1/2) due to the symmetry. Thus, in the following, 
we will consider Na = Nb- Above the critical density p c , 
the system separates in two phases I and II. One phase 
is rich in A particles and the other is rich in B parti- 
cles. Furthermore, due to the symmetry, those phases 
are symmetric in composition (x A = x B and x A = x B ). 

The determination of the coexistence curve x A and x l A l 
is possible from the use of the two equivalent boxes. The 
overall composition xa — Na/{Na + Nb) is fixed dur- 
ing the simulations whereas the particular compositions 
x A = N{/(N A + N B ) and x l A l = N 1 / / (N 1 / + N B r ) of the 
boxes I and 77 are free to fluctuate. N^ corresponds to 
the number of particles of component X = A or B inside 
the box Y — I or II. Notice that the number of particles 
in each box N 1 = N A + N B and N 11 = N A + N B also 
fluctuates. The ensemble considered is thus intermediate 
between the grand canonical ensemble and the semigrand 
canonical ensemble. We will discuss in the following the 
consequences of the slight density fluctuations inside each 
box. 

In complement to the number density p, it is interest- 
ing to introduce the scaled packing fraction 7/ = vabP 
with vab the volume of unlike particles, vab — 7T o- AB /4 



in two dimensions and ira AB /6 in three dimensions. This 
definition of the packing fraction which takes into account 
the non-additivity allows us to compare the critical pack- 
ing fraction of the WR model with those of the NAHC 
mixtures. Notice that from this definition, the packing 
fraction may exceed one. 



B. Description of the cluster algorithm 

A cluster algorithm has been introduced recently by 
Dress and Krauth |18| for the simulation of hard-core 
mixtures. Inspired by the lattice cluster algorithms of 
Swendsen-Wang [2^] and Wolff 23], it allows one to 
equilibrate large off-lattice systems (up to 10 6 particles) 
thanks to a non-local move of a large number of particles 
at each Monte-Carlo step. The general idea of the algo- 
rithm is to take advantage of the hard-core interactions 
between particles to construct a cluster of particles that 
will be moved at each Monte-Carlo step satisfying the 
detailed balance while keeping ergodicity. 

The Monte-Carlo step is constructed as follows: 

i) We select randomly one of the two boxes. 

ii) A second box is chosen randomly (it can be the 
same as the previous one). 

iii) An inversion symmetry around a randomly chosen 
pivot point is performed on all the particles of the second 
box. 

iv) The two boxes with their particles are then super- 
imposed on top of each other resulting in a set of clusters 
of overlapping particles (in the sense of the hard-core in- 
teractions) . 

v) A particle is randomly selected and the cluster from 
which it belongs is flipped. Each particle belonging to 
this cluster is moved from its initial box to the other 
box in the position corresponding to the superimposed 
configuration. 

For a clear graphical representation of a Monte-Carlo 
step see Figure 1 in Dress and Krauth paper 0] . 

Let us first prove that the new configuration obtained 
after the Monte-Carlo step satisfies all the hard-core in- 
teractions between particles, a) For two particles outside 
the flipped cluster there is no overlap since those parti- 
cles did not move and there was no overlap before the 
Monte-Carlo step, b) For two particles inside the clus- 
ter, they were flipped keeping there relative positions in 
such a way that there is still no overlap, c) Finally, by 
construction of the cluster, there is no overlap between a 
particle belonging to the cluster and one outside the clus- 
ter when the two boxes are superimposed and a fortiori 
after the Monte-Carlo step. Those arguments justify the 
fact that the new configuration obtained satisfies all the 
hard-core interactions between particles in each box. 

Now let us justify the detailed balance of the cluster al- 
gorithm. Each Monte-Carlo step has its symmetric step 
in the sense that if we start from the final configuration 
there is a pivot point that gives back the initial config- 
uration. The pivot points being randomly selected with 
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a uniform distribution, there are equality between the 
probabilities to go from the initial to the final configura- 
tion and vice- versa. The detailed balance is thus satisfied 
noting that due to the hard-core interactions both con- 
figurations have the same equilibrium probability. In the 
case of a general potential of interaction between parti- 
cles, the cluster algorithm may be generalized with the 
extent of some rejections of the Monte-Carlo steps |24j . 

On the point ii) of a Monte-Carlo step, either the same 
box or the other one may be selected. This leads to 
intra-box or inter-box Monte-Carlo steps. The reason 
for this choice is to decrease the equilibration time and 
as we will see to satisfy the ergodicity of the cluster al- 
gorithm. With intra-box Monte-Carlo steps (same box 
selected twice), it has been shown already 113 that 
the algorithm satisfies internal ergodicity like the usual 
Monte-Carlo algorithm. By internal ergodicity, we con- 
sider the fact that for the given composition of the box, 
all possible configurations of the particles are attainable 
by the cluster algorithm. In fact, if the pivot point is cho- 
sen sufficiently close to a particle position and this par- 
ticle considered as the starting point for the cluster con- 
struction, the Monte-Carlo step corresponds to a slight 
move of this particle without affecting the other particles. 
This move corresponds to a usual Monte-Carlo move in 
a general algorithm. This argument justifies the internal 
ergodicity. However, from those moves, the composition 
of the box is kept constant. In order to change the compo- 
sition or the relative number of particles of components A 
or B in each box, inter-box Monte-Carlo steps are neces- 
sary. Those inter-box Monte-Carlo steps are performed 
to exchange the components of particles. If the pivot 
point is selected such that two particles from different 
components and boxes superimposed exactly and if one 
of those particles is selected as the cluster starting point, 
the Monte-Carlo step reduces to the exchange of those 
particles. The Monte-Carlo move then corresponds to 
the usual change of components considered by the Monte- 
Carlo algorithm in the semigrand canonical ensemble. In 
summary, both usual Monte-Carlo steps (the move of 
a particle and the change of components) are possible 
moves in this cluster algorithm justifying the ergodicity 
if this ergodicity is assumed for the usual Monte-Carlo 
algorithm in the semigrand canonical ensemble. 

As can be seen from the discussion on the ergodicity, 
the use of the two boxes is useful for the non-additive 
hard-core mixtures. It has another strong advantage 
since it allows one to determine the coexistence curve 
or the relative composition of the two separated phases 
above the critical density. Due to the symmetry of the 
problem, the critical point corresponds to an equal parti- 
tion in particles A and B {x c — 1/2) and above the crit- 
ical density the coexistence curve is symmetric around 
x c = 1/2. The use of identical boxes is thus justified 
since the densities p 1 and p 11 of the two phases are equal. 
As previously discussed, the total number of particles in 
each box is slightly fluctuating with this cluster algo- 
rithm in contrary to the simulations in the semigrand 



Exponent 


V 


7 


P 


2D (exact) 


1 


7/4 


1/8 


3D (num.) 


0.627 


1.239 


0.326 



TABLE I: Critical exponents for the Ising universality class 
in two and three dimensions. 



canonical ensemble. However, those fluctuations of the 
density are not related to the composition fluctuations in- 
side each box. More importantly, around the phase sep- 
aration transition, the composition fluctuations diverge 
whereas the density fluctuations stay insensitive to the 
transition. As a consequence the slight density fluctua- 
tions do not affect the coexistence curve determined from 
the cluster algorithm. 

The last question concerning the cluster algorithm con- 
cerns the equilibration time or the number of Monte- 
Carlo steps necessary to equilibrate the system. The 
Swendson-Wang cluster algorithm was introduced to sim- 
ulate systems of Ising spins with ferromagnetic interac- 
tions between neighbor spins on a lattice [22| . In contrary 
to the simple Monte-Carlo algorithm, this cluster algo- 
rithm does not suffer from a critical slowing down at the 
phase transition due to the fact that the flipped clusters 
are then directly related to the spin clusters observed 
around the phase transition. In the case of the present 
cluster algorithm, such direct relation of the flipped clus- 
ters with the configurations of particles is not demon- 
strated. However, the number of Monte-Carlo steps nec- 
essary for the equilibration of the system does not seem 
to increase significantly when approaching the phase sep- 
aration transition. It is also interesting to note that this 
number is roughly independent of the system size as ob- 
served for cluster algorithms on lattices and in contrary 
to usual Monte-Carlo algorithms where this number usu- 
ally increase strongly with the system size. Another im- 
portant point is that there is no critical slowing down for 
the equilibration when the non-additivity is increased. 
Thus, this cluster algorithm is well adapted for large non- 
additivities A and for the WR model in comparison to 
the Monte Carlo simulations in the semigrand canonical 
ensemble. In fact, the critical slowing down is observed 
for small non-addtivities when the phase separation tran- 
sition occurs at high densities close to a fluid-solid tran- 
sition. 



III. RESULTS OF THE NUMERICAL 
SIMULATIONS 



A. Critical exponents and critical packing fractions 

As all second order phase transitions, the critical point 
in NAHC mixtures is characterized by different critical 
exponents. From finite size scaling of equilibrium ther- 
modynamic quantities close to the critical point (p c ,x c ), 
it is possible to determine those exponents. For the 
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NAHC mixtures and for the WR model, the critical ex- 
ponents are supposed to belong to the Ising universality 
class (see numerical values of v, 7 and (3 in Table . In 
the following, we define different ways to determine or 
test the value of the critical exponents. We also describe 
four different ways to define a finite size critical packing 
fraction r} c (L) from numerical simulations. 

From finite size scaling, it is possible to extract the 
(infinite size) critical packing fraction rj c [T2I |25| : 



T) C (L) = rj c - A/L 



l/u 



(1) 



The finite size of the system L is defined as Na = L 
where d is the dimension of the system. The critical ex- 
ponent v is independent of the definition of the finite size 
critical packing fraction considered in contrary to the co- 
efficient A. Due to the corrections to scaling for small 
system sizes, it is usually difficult to extract the criti- 
cal exponent v. However, a linear behavior of the finite 
size critical packing fraction with the rescaled system size 
1/L 1 / 1 ' is a good test of the universality class of the model 
considered and allows one to extract the critical packing 
fraction r\ c . 

The phase separation transition is characterised by the 
order parameter m — 2{xa — x c ) with — 1 < m < 1. Its 
probability distribution P(m; 77, L) for a packing fraction 
77 and a system size L changes form around the critical 
packing fraction rj c . In the thermodynamic limit (L = 
00), the distribution of the order parameter presents 
delta picks. It has a single pick at m = below the 
critical packing fraction (77 < rj c ). However, above the 
critical packing fraction (77 > rj c ), the distribution is dou- 
ble picked at values ±T7i max with m max ~ (77 — rj c ) 13 for 
V ^ Vc- For finite size systems, the picks broaden but 
the change of the distribution remains and a maximum 
(77, L) > may be defined for each packing fraction 
77 and size L above the size dependent critical packing 
fraction 77™ ax (L). The dependence of the maximum m max 
on the packing fraction at a fixed system size is pfij : 



m max ~ ( V - V^Wf 



for -q > r]™ ax (L). Knowing the exponent /3, it is thus 
possible to extract the finite size critical packing fraction 
i)™(L) from a linear fit of ro^x as function of r/. For 
high packing fraction 77, m max saturates to one. This 
limits the range of packing fractions for which the lin- 
ear fit is valid. Due to this limited range of the power 
law behavior, it is difficult to extract the critical expo- 
nent (3. However, a linear dependence obtained for the 
expected value of the critical exponent (3 is still a strong 
confirmation of the universality class. 

Due to the symmetry of the model, the average of the 
order parameter is zero. However, the mean absolute 
value of the order parameter is another possibility to de- 
fine a finite size order parameter: 



(\m\)(r,,L) 



\m\P(m; 77, L)dm. 



(3) 



A 


0.5 


1.0 


2.0 


00 


Ising 


2D 


1.749(8) 


1.749(8) 


1.746(7) 


1.742(8) 


7/4 


3D 


2.03(8) 


2.05(8) 


2.05(8) 


2.05(8) 


1.98 



TABLE II: Numerical results for the ratio of critical expo- 
nents 7/1/ for different non-additivities and for the WR model 
(A = 00) in two and three dimensions. Numbers in paren- 
thesis correspond to the error on the last digits. The data 
in column Ising are the predictions of the Ising universality 
class. 



This definition suffers from an additional drawback com- 
pared to m max . The average absolute value of the order 
parameter does not vanish at the finite size critical pack- 
ing fraction 77™ (L). Thus, the power law behavior [2o| : 



<|m|)^(77 - 77r(i)) /3 



(4) 



presents corrections to scaling not only for large packing 
fraction but also around 77™ (L). This finite size criti- 
cal packing fraction 77™ (L) may still be extracted from a 
linear fit of ((ml) 1 /' 3 on a limited range of packing frac- 
tions. The small value of the critical exponent (3 in two 
dimensions leads to a sufficiently large range, however, in 
three dimensions, the larger value of (3 renders the range 
of power law behavior Q too small to be able to extract 
V ™(L). 

The maximum of the modified susceptibility ^(77, L) = 
(m 2 ) — (|tti|) 2 is a third possibility to define a finite size 
critical packing fraction. This modified susceptibility 
presents a single maximum in contrary to the real sus- 
ceptibility (m 2 ) — (to) 2 . The general form of the modified 
susceptibility is |25j : 



X (v,L)^L''^- d x(L 1 ^(r 1 -r, c )) 



(5) 



with x(x) a function with a single maximum Xmax- The 
packing fraction at the maximum of the modified sus- 
ceptibility defines the finite size critical packing fraction 



(2) Vc(L). The modified susceptibility may also be used to 
extract the ratio of critical exponents 7/z/. Its maximum 
depends algebraically on the system size with an expo- 
nent 7/1/ -d [25|: 



Xmax (77 c *(L),L)~^/- 



d Y 
A max- 



(6) 



The ratio of critical exponents 7/1/ is thus simply deter- 
mined by the slope of a linear fit in a log-log scale. The 
results for the two and three dimensional systems are 
presented in Table [H] in comparison with the Ising uni- 
versality prediction. Those results are discussed later. 

Another possibility to determine the critical packing 
fraction concerns the Binder parameter [25|: 



U(r),L) = l- 



3(m 2 ) 



(7) 



The Binder parameter saturates to 2/3 for large packing 
fractions and vanishes for small ones. It also presents the 



interesting property to intersect at the critical packing 
fraction at least for sufficiently large system sizes. The 
value U(r] c , L) = U* at this intersection is expected to be 
universal. The monotonous behavior of the Binder pa- 
rameter may be used to define a finite size critical packing 
fraction r]f(L) as the value for which U(r]^(L), L) = 1/2. 
The choice of the value 1/2 is arbitrary and could be mod- 
ified as soon as it is sufficiently different from the bound- 
ary values and 2/3 and from the intersction value U* . 
The particular choice of U* would render the finite size 
packing fraction independent of the system sizes. This 
point will be discussed further for the three dimensional 
systems. 

We defined different critical packing fractions for finite 
size systems from the maximum of the modified suscepti- 
bility ri%(L), from the Binder parameter rj^(L) and from 
the coexistence curve either from the maximum of the 
distribution of the order parameter 7y™ ax (L) or from the 
average absolute value of the order parameter r/^ v (L). 
The numerical results obtained from the cluster algo- 
rithm are analyzed in the following sections. 



B. Two dimensional NAHC mixtures 

In the two dimensional model, we considered four dif- 
ferent non-additivities A = 0.5, 1.0, 2.0 and 4.0 as well as 
the WR model (A = oo). The system sizes ranged from 

1/2 

L = = 10 to 400. The largest systems contained 

Na+Nb — 320000 particles. Numerical simulations were 
divided in five consecutive runs of 10 5 to 10 6 Monte-Carlo 
steps depending on the system sizes with longest runs for 
the largest systems. The first run is kept for equilibration 
of the initial configuration and the last four runs for the 
data collection and error estimation. 

On Fig. ^ we plot equilibrium configurations for two 
different packing fractions. The configurations corre- 
spond to systems with Na + Nb = 3200 particles and 
a non-additivity A = 2. The likewise hard-core diame- 
ter a aa and obb are represented respectively by white 
and black disks whereas the unlike diameter a ab is rep- 
resented by light and dark grey disks respectively for the 
A and B particles. No overlaps between unlike particles 
are present but overlaps of likewise particles are observed 
on the grey scale. However, no overlap is present be- 
tween black and white particles justifying that both con- 
figurations satisfy all the hard-core interactions. For the 
top configuration above the phase separation transtion, 
i] = 1.03 > T] C (L) ~ 0.95, we observe a large difference 
between the number of A and B particles. The large per- 
colating cluster of B particles (the dark grey particles) 
is a clear evidence of the phase separation. On the con- 
trary, on the bottom configuration, r\ = 0.72 < r) c (L), 
the A and B particles are roughly in equal number and 
perfectly mixed. In this systems, the phase separation 
did not occur, however, due to the stronger unlike parti- 
cles hard-core interactions, a local clustering of likewise 
particles is present. 
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FIG. 1: Equilibrium configurations for a system above (top) 
and below (bottom) the phase separation transition. The 
systems correspond to Na + Nb ~ 3200 particles with a non- 
additivity A = 2. White and black disks correspond to like- 
wise diameter cjaa and <jbb whereas the light and dark grey 
disks correspond to the unlike diameter gab- With this rep- 
resentation, the overlaps are only allowed by the hard-core 
interactions to the likewise particles on the grey scale disks. 



The modified susceptibility x is determined as func- 
tion of the packing fraction rj for different system sizes. 
It presents a single maximum XmaxiVci-L), L) at the fi- 
nite size critical packing fraction rj^(L). As can be seen 
on Fig. |3 this maximum depends on the system size with 
a power law behavior. From a linear fit in a log-log scale, 
we deduce the ratio of critical exponents j/p (see Ta- 
ble Due to the corrections to scaling observed for 
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FIG. 2: Maximum of the modified susceptibility Xmax as func- 
tion of the system size L in a log-log scale for three different 
non-additivities A and for the WR model ( Widom) . The ratio 
of critical exponents 7/1/ is determined from linear fits. The 
three smallest system sizes are removed from the linear fit due 
to the corrections to scaling effects for those small sizes. The 
data have been shifted for clarity. 



small system sizes, the three smaller sizes (L — 10, 12 
and 15) where removed for the linear fit to extract 7/f. 
This ratio of critical exponents compare nicely with the 
Ising universality class prediction 7/4 for all the non- 
additivities considered and for the WR model. The esti- 
mation for the errors presented on Tabic ITT1 comprises the 
error on the maximum of the modified susceptibility and 
the error coming from the linear fit. The relative error 
on the ratio of critical exponents 7/1/ is smaller than 1%. 

The confirmation of the critical exponent is obtained 
from the two different definitions of the finite size or- 
der parameter: m max and (|m|). First, we plot the 
rescaled maximum of the distribution of the order pa- 
rameter mmax as function of 77 for different system sizes 
and for a non-additivity A = 2 on Fig. Second, we 
plot the rescaled average order parameter (|m|) 1 /' 3 for 
the same system sizes but for a non-additivity A = 1 on 
Fig.yp. The linear behavior observed on both figures for 
a critical exponent 0=1/8 confirms the Ising universal- 
ity class prediction. Similar results are obtained for all 
the non-additivities and for the WR model. The range 
of the linear regime is still rather small and limited for 
large packing fractions to mUaL < 0.7 and (|m|) 1 /' 3 < 0.6. 
Furthermore, due to the strictly positive average of the 
absolute order parameter, the power law behavior is also 
limited from below to (|m|) 1 /' 3 > 0.1. The restriction 
on the range of the power law behavior prevents from a 
direct determination of the critical exponent 0. 

The critical packing fraction r\ c is then determined 
from the plot of the finite size critical packing fraction 
r] c (L) as function of the rescaled system size 1/ 'L x l v with 
the expected exponent v = 1. On Fig. the numer- 




0.8 0.85 0.9 0.95 

FIG. 3: (a) Rescaled maximum of the distribution of the 
order parameter rnUaL for a non-additivity A = 2 and (b) 
rescaled average order parameter (|m|) 1///3 for a non-additivity 
A = 1 as function of the packing fraction rj for different system 
sizes. Linear fits confirm the critical exponent 0=1/8 pre- 
dicted by the Ising universality class and allow us to extract 
the finite size critical packing fractions ??™ ax (L) and r]c V (L) 
respectively. 



ical results for ry™ ax , 77™ and rj^ are plotted for a non- 
additivity A = 1. The results for 77* close to those for rj^ 
have been removed for clarity. The error bars are smaller 
than the symbols and thus not represented on the figure. 
As can be seen on Fig. 0] ?y™ ax and 77°" present a linear 
behavior with respect to 1/L 1 ^ for v = 1. This confirms 
the value of the critical exponent v predicted by the Ising 
universality class. However, rj^(L) presents corrections 
to scaling and a quadratic fit is necessary to extract the 
critical packing fraction r)f? . The same is true for 77* and 
similar results are obtained for all the non-additivities as 
well as the WR model. 

All the critical packing fractions from the four different 
definitions of their finite size analog arc presented on Ta- 
ble [ffl] for all the non-additivities considered and for the 
WR model. It is interesting to notice that the four differ- 
ent definitions are not identical since the numerical val- 
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Their critical packing fraction obtained for A = 1 is 
7y c = 0.118(1 + A) 2 = 0.472. A scaled particle theory 
predicts rj c = 0.096(1 + A) 2 = 0.385 H3 whereas a virial 
expansion 0.324. The different definitions of the packing 
fraction in both cases explains the A dependence intro- 
duced to compare to our prediction r\ c = 0.897(2). The 
strong difference should be contrasted by the fact that 
the prediction of a first-order perturbation theory usu- 
ally referred to as MIX1 gives higher values for the criti- 
cal packing fractions especially for large non-additivities 
(A > 0.4) [l(j]. It could be interesting to pursue the com- 
parison of this theory for higher non-additivities (A > 1) 
with our numerical predictions. 



C. Three dimensional NAHC mixtures 



FIG. 4: Finite size critical packing fractions rj^ v (L), ^^(L) 
and r]c (L) as function of the rescaled system size 
The non-additivity considered is A = 1. Linear fits of rf^ v 
and 77™ ax confirm the critical exponent v = 1 of the Ising 
universality class and allow us to extract the critical packing 
fractions n^f° and 7/™ ax . A quadratic fit is necessary for rj^{L) 
in order to extract the critical packing fraction 77^ due to the 
corrections to scaling and stronger finite size effects. 



A 


c 


max 
'lc 






0.5 


0.8141(5) 


0.8138(11) 


0.811(2) 


0.811(4) 


1.0 


0.8988(8) 


0.8976(20) 


0.896(4) 


0.895(5) 


2.0 


1.0215(18) 


1.0179(45) 


1.017(9) 


1.017(6) 


4.0 


1.147(4) 


1.141(5) 


1.141(5) 


1.140(5) 


CO 


1.228(3) 


1.222(6) 


1.224(8) 


1.225(5) 



TABLE III: Numerical results for the critical packing fraction 
for the two dimensional NAHC mixtures for different non- 
additivities and for the WR model (A = 00). The numbers 
in parenthesis correspond to the error on the last digits. 



ues obtained for finite sizes differ (see Fig. 0J. However, 
the thermodynamic limit results for the critical packing 
fractions are consistent for all the four definitions. The 
estimation of the relative error which is smaller than 1% 
for all the critical packing fractions combines the error 
on the finite size estimates and the error coming from 
the finite size scaling to extract the thermodynamic limit 
results. 

Let us now compare our results to previous simula- 
tions. Our prediction for the critical density Pc&\b = 
1.560(10) of the WR model is in good agreement with 
the value 1.566(3) obtained by Johnson et al. |3j. For the 
NAHC mixtures in two dimensions, the determination of 
the critical packing fractions have been recently obtained 
from numerical simulations by Saija and Giaquinta [To| . 
Only small non-addtivities A < 1 have been considered 
and the system was limited to 800 particles. The critical 
packing fraction was determined without finite size scal- 
ing and may thus be considered only as a lower bound. 



In three dimensions, we considered a large number of 
different non-additivities A = 0.5, 0.6, 0.8, 1.0, 2.0, 3.0, 

4.0 and 9.0 as well as the WR model (A = 00). The 

1/3 

system sizes ranged from L = N A = 5 to 50. The 
largest systems contained Na + Nb = 250000 particles. 
Systems with up to 2 x 10 6 (or L = 100) where simulated 
for some non-additivities. The equilibrium time was then 
close to the day on simple personal computers. Since a 
large number of densities needs to be simulated to extract 
the critical density, we are for L — 100 at the limit of 
the numerical capabilities. As for the two dimensional 
systems, five consecutive runs were simulated. The first 
one was used to equilibrate the initial configuration and 
the four remaining ones served for the collection of data 
and the estimation of errors. 10 5 to 10 6 Monte-Carlo 
steps for each run were sufficient for the equilibration 
and in order to obtain a small error. 

The modified susceptibility x is determined as function 
of the packing fraction r\ for different system sizes. It 
presents a single maximum x max (7y*(L), L) at the finite 
size critical packing fraction rj^(L). As can be seen on 
Fig. [31 this maximum depends on the system size with a 
power law behavior. From a linear fit of the maximum of 
the modified susceptibility as function of the system size 
in a log- log scale, we deduce the ratio of critical exponents 
j/v (see Table Hl*|) . Those ratios for all non-additivities 
are systematically larger than the Ising universality class 
prediction but still inside the error bars. The relative 
errors were estimated as the sum of the average relative 
error on the maximum x max and the error due to the 
linear fit. 

The existence of corrections to scaling may explain 
the over-estimation of "f/v. Removing the data from 
the small size systems for the linear fits reduces the dis- 
crepancy with the Ising universality class prediction. A 
possible source for the corrections to scaling observed is 
the existence of fluctuations on the density inside each 
box during the simulations. Even though those fluctu- 
ations are small they may affect slightly the scaling be- 
havior especially for small systems were the fluctuations 
are stronger. Some constrained simulations have been 
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FIG. 5: Maximum of the modified susceptibility x 

max £LS 

function of the system size L in a log-log scale for different 
non-additivities A and for the WR model ( Widom) . The ratio 
of critical exponents ■y/v is extracted from linear fits. Data 
have been shifted for clarity. 
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FIG. 7: Finite size critical packing fractions 7/Jf ax (L) (a) 
and r/*(L) (b) as function of the rescaled system size 
with v — 0.627 as predicted by the Ising universality class. 
Different non-additivities are considered as well as the WR 
model (same legends for both figures). Linear fits of ?7™ ax as 
well as rj* confirm the critical exponent v = 0.627 and allow 
us to extract the critical packing fractions 77" lax and rf£ . 



FIG. 6: Rescaled maximum of the distribution of the order 
parameter mi/ax as function of the packing fraction r\. The 
non-additivity considered is A = 0.8. The critical exponent 
13 — 0.326 predicted for the Ising university class has been 
used and the critical packing fractions ?7™ ax (L) for different 
system sizes L are extracted from linear fits. 



done were the density in each box was kept constant. 
The clusters were moved only if they did not change the 
total number of particles inside each box. This modifi- 
cation of the algorithm leads to similar results with less 
pronounced corrections to scaling. However, errors were 
also stronger since due to the rejections the equilibration 
time was sensibly increased. 

The confirmation of the critical exponent fj is obtained 



from the plot of the rescaled maximum of the distribution 
of the order parameter ml^ as function of the packing 
fraction r\ for different system sizes. A non-additivity 
A = 0.8 has been considered on Fig. but similar re- 
sults are obtained for all non-additivities and for the WR 
model. The linear behavior observed on the figure for a 
critical exponent (3 — 0.326 confirms the Ising univer- 
sality class prediction. The range of the linear regime is 
still rather small and limited for large packing fractions to 

m\LL < 0.5. The range of the linear regime for (jml) 1 '^ 
is even smaller and cannot allow us the confirmation of 
the critical exponent (3 neither the determination of the 
finite size critical packing fraction rf^"{L). 

The critical packing fraction r\ c is determined from the 
plot of the finite size critical packing fraction T) C (L) as 
function of the rescaled system size 1/L 1 /" with the ex- 
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A 


max 


rg 




0.5 


0.3577(3) 


0.3574(3) 


0.3574(3) 


0.6 


0.3572(4) 


0.3564(4) 


0.3568(4) 


0.8 


0.3583(5) 


0.3581(4) 


0.3580(5) 


1.0 


0.3610(5) 


0.3614(3) 


0.3614(5) 


1.5 


0.3706(5) 


0.3703(5) 


0.3705(8) 


2.0 


0.3758(5) 


0.3766(5) 


0.3766(8) 


3.0 


0.3827(10) 


0.3839(6) 


0.3843(9) 


4.0 


0.3862(10) 


0.3871(8) 


0.3872(11) 


9.0 


0.3896(9) 


0.3908(6) 


0.3914(9) 


oo 


0.3910(4) 


0.3912(4) 


0.3912(4) 



TABLE IV: Numerical results for the critical packing frac- 
tions in three dimensions for the different non-additivities con- 
sidered and for the WR model (A = oo). The numbers in 
parenthesis correspond to the error on the last digits. 



pected exponent v — 0.627. The numerical results for 
i) c max (L) (Fig. Ek) and for rj*{L) (Fig. Eb) are plotted 
for different non-additivities and for the WR model. The 
error bars are smaller than the symbols and thus not rep- 
resented. As can be seen, ^™ ax and rf£ present a linear 
behavior with respect to 1/ ' L x l v for v = 0.627. This con- 
firms the value of the critical exponent v predicted by 
the Ising universality class. The results are similar for 
rjc and for the other non-additivities considered. The 
corrections to scaling observed for 77* and r/^ in two di- 
mensions are not present in the three dimensional sys- 
tems. It is possible to notice a change in the slope of 
the linear fit around A = 1. Already evident on Fig. [7] 
this is confirmed by the results for the non-additivities 
A = 0.6 and 0.8 (data not shown on the figure). For 
A > 1 the slope is independent of the non-additivity 
whereas for A < 1 the slope is increasing with the non- 
additivity. This change may be due to the onset of the 
fluid-solid transition that the model experience for large 
packing fractions. At sufficiently low non-additivity, the 
phase separation transition is expected to be pre-empted 
by the fluid-solid transition J27|. 

All the critical packing fractions for the different non- 
additivities and for the WR model are presented on Ta- 
ble El The relative errors are less than 0.5% and the 
three definitions of the critical packing fraction ry™ ax , tj * 
and rjc lead to identical estimations in the thermody- 
namic limit. 

Let us now compare our results with others. Our 
prediction for the critical density of the WR model is 
Pc a AB = 0.7470(8) which is in strong agreement with 
the value 0.748(2) obtained by Johnson et al. |3j with an- 
other cluster algorithm. Our result improves by a factor 
two the error on this critical density. The critical pack- 
ing fraction has been determined from different numerical 
simulations for a large set of non-additivities. For a non- 
additivity A = 1, our prediction ry c /(l+A) 3 = 0.04515(8) 
is slightly higher than the value 0.04484(20) proposed 
by Gozdz |l2j. However, it strongly confirms the previ- 
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FIG. 8: Critical packing fraction for the different non- 
addtivities r? c (A) as function of 1/(1 + A) 2 (disks) and for 
the WR model (square). The line corresponds to a linear fit 
for the six largest non-additivities and leads to an estimation 
of rjc = 0.3907(5) for the WR model. 



ous under-estimation of 0.0288 by Saija et al. |9j. The 
non-additivity dependence introduced is due to a dif- 
ferent definition for the packing fraction. Gozdz [T^ 
also predicted the critical packing fraction 0.0866(2) and 
0.0611(2) for the non-additivity A = 0.6 and 0.8 respec- 
tively. Its predictions compare nicely with our results 
r? c /(l + A) 3 = 0.0871(1) and 0.0614(1). It is important 
to notice that the critical packing fraction is determined 
to our knowledge for non-additivities A > 1 for the first 
time in this study. This is possible thanks to the absence 
of a critical slowing down of the simulations for large non- 
additivities with the cluster algorithm in comparison to 
the simulations in the semigrand canonical ensemble. 

From the determination of the critical packing fraction 
for large non-additivities, it is possible to consider the 
convergence of the NAHC mixtures to the WR model 
which corresponds to an infinite non-additivity. It is 
evident from Fig. [S] that the critical packing fraction 
rj c (A) = i] c (oo) — B / (1 + A) 2 for large non-additivities. 
A linear fit for the six largest non-additivities leads to 
an estimate of r\ c — 0.3907(5) for the WR model in good 
agreement with the direct simulations. The coefficient B 
was estimated to be 1.21. This interesting dependence of 
rj c (A) on A may be a good test for the different theories 
developed to determine the critical packing fraction of 
the NAHC mixtures. 

The Binder parameter offers another test of the uni- 
versality class from its intersection value U* . In Fig. 
we have plotted the Binder parameter as function of the 
packing fraction for different system sizes and for the 
non-additivity A = 2. It can be seen that the intersec- 
tions between the different curves occur at smaller and 
smaller values as the system size increases but still at 
higher values than the expected universal one U* = 0.47. 
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FIG. 9: Binder parameter for different system sizes as func- 
tion of the packing fraction r\ for the non-additivity A = 2.0. 
The horizontal line represents the universal intersection value 
U* = 0.47 expected for the Ising universality class. 



Similar results are observed for the other non-additivities 
and for the WR model. This result was already observed 
in Gozdz simulations |l2j |. Corrections to scaling may 
explain such behavior. For our definition of the critical 
packing fraction r)J? based on the Binder parameter, the 
arbitrary choice was U(r)f(L),L) = 1/2 > U* . Thus, 
rjf(L) was expected to decrease with increasing system 
size. The opposite behavior is found. The absence of 
corrections to scaling in the dependence of r]®(L) with 
the system size is thus an argument against this expla- 
nation for the high intersection in Fig. El However, no 
other explanations have been found. 



for the WR model. This cluster algorithm allows one 
to simulate large systems (up to 10 6 particles). Each 
Monte-Carlo step corresponds to the non-local move of a 
large number of particles reducing strongly the equilibra- 
tion time. The absence of critical slowing down for in- 
creasing non-additivities A allows us to study large non- 
additivities and the convergence of the NAHC mixtures 
to the WR model when A — > oo. 

Two and three dimensional systems have been consid- 
ered. In both cases, the models were found to belong 
to the Ising universality class for all the non-additivities 
considered as well as for the WR model. The critical 
packing fractions r\ c have been determined from four dif- 
ferent finite size definitions. All the definitions lead to 
identical predictions in the thermodynamical limit with 
a high precision. 

Different theories have been used to determine the 
critical packing fractions T) C (A) of NAHC mixtures [26| . 
from the scaled particle theory [2]} to the virial expan- 
sion [2i| or the density functional theory [^3]. A 
comparison with our precise numerical predictions of the 
critical packing fractions for a large set of non-additivities 
A is possible. Even if a quantitative comparison between 
theory and simulations is difficult, the A behaviour of 
the critical packing fraction (j] c = ?y c (co) — BJ (1 + A) 2 in 
three dimensions) for large non-additivities A is a good 
test for a theory. 

Recently, Jagannathan and Yethiraj 5] have studied 
the dynamical behavior of the WR model close to the 
phase separation transition. The cluster algorithm could 
be used to consider larger systems in order to be closer 
to the critical density. The cluster algorithm could also 
be used to study the stability and interfacial properties 
in confined geometries [3l|. 



IV. DISCUSSION 
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